#ifndef SYSUTILS_H
#define SYSUTILS_H


class SysUtils
{
 public:
  static TChain * InputFiles(const TString file, const TString tr, Char_t *dir=0x0);
  static void ToNaturalScale(TAxis *ax);
  static void ToNaturalScale(TH1 *hh);
  static void ToNaturalScale(TGraphAsymmErrors* gr);

  static Int_t GetCommonLimits(const Int_t ngr,  const TGraphAsymmErrors * grs[], Int_t id0s[], Int_t id1s[]);
  static TGraphAsymmErrors * HistToGraph(const TString grname, const TH1 *hh, const Double_t thres=0, const TH1 *herr=0x0, const Double_t xmin=-1e10, const Double_t xmax=1e10, const Bool_t klogx=kFALSE);

  static TGraphAsymmErrors * TreeToGraph(const TString tname, const TString cut, Char_t * sx, Char_t * sy, Char_t * sex0=(Char_t*)"", Char_t * sex1=(Char_t*)"", Char_t * sey0=(Char_t*)"", Char_t * sey1=(Char_t*)"", const Bool_t kXLogToLinear=kFALSE, const Bool_t kprint=kFALSE);

  static void FitSlicesY(const TH2D *hh, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, TH1D *&hchi, const Double_t thres, const Double_t lowfrac, const Double_t highfrac, TList *ll);
  static Double_t ToPt(const Double_t mom, const Double_t eta);
  static void ScaleGraphPoint(TGraphAsymmErrors * gr, const Int_t ipoint, const Double_t fac);
  static TGraphAsymmErrors * MeasurementConversion(const TString hn, const TGraphAsymmErrors * gv1, const TH2D * hy2x1, const Bool_t kcutEdge, const Double_t thres);

  static void GetGrSysMode0(TList *ll, const TString tag, const Bool_t *badgr, const Bool_t kfit, const Double_t hmin, const Double_t hmax);
  static void BinLogX(TAxis *axis);

 private:
  static Bool_t GetGrPointLimits(const TGraphAsymmErrors * gr, Double_t & gmin, Double_t & gmax, Bool_t & klogx);
  static Double_t * GetAxisArray(TAxis * aa);
  static Double_t TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr);
  static Double_t GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms, Double_t * gerr);
  static Double_t TruncatedMean(const Int_t nx, const Double_t xdata[], const Double_t lowfrac, const Double_t highfrac, Double_t * grms=0x0, Double_t * gerr=0x0, Double_t *wws=0x0);

};

TChain * SysUtils::InputFiles(const TString file, const TString tr, Char_t *dir)
{
  TChain *ch=new TChain(tr);

  if(file.Contains(".root"))
    ch->Add(file);
  else{
    ifstream fin(file);
    if(!fin){
      printf("SysUtils::InputFiles file not found \n%s\n\n",file.Data()); exit(1);
    }

    TString buff;
    while(fin.good()){
      fin>>buff;
      if(buff!=""){
        if(dir){
          buff.Prepend(dir);
        }
        ch->Add(buff);
      }
    }
  }

  //const Int_t ent = ch->GetEntries(); //takes infinity time!!                                                         
  printf("\t%d trees!\n",ch->GetNtrees());

  return ch;
}


void SysUtils::ToNaturalScale(TAxis *ax)
{
  TAxis* oldx = (TAxis*)ax->Clone("oldx");
  ax->SetLimits(TMath::Power(10,oldx->GetXmin()), TMath::Power(10,oldx->GetXmax()));
  const Int_t nb = oldx->GetNbins();
  Double_t *bins = new Double_t[nb+1];
  bins[0]=TMath::Power(10,oldx->GetXmin());
  for(Int_t ii=1; ii<=nb; ii++){
    bins[ii]=TMath::Power(10,oldx->GetBinUpEdge(ii));
  }
  ax->Set(nb, bins);

  delete oldx;
  delete bins;
}

void SysUtils::ToNaturalScale(TH1 *hh)
{
  ToNaturalScale(hh->GetXaxis());
}


void SysUtils::ToNaturalScale(TGraphAsymmErrors* gr)
{
  for(Int_t ip=0; ip<gr->GetN(); ip++){
    const Double_t xx = gr->GetX()[ip];
    const Double_t exl = gr->GetEXlow()[ip];
    const Double_t exh = gr->GetEXhigh()[ip];

    const Double_t yy = gr->GetY()[ip];
    const Double_t eyl = gr->GetEYlow()[ip];
    const Double_t eyh = gr->GetEYhigh()[ip];

    //printf("%d %e %e\n",ip,xx,yy);                                                                                                                                                                
    gr->SetPoint(ip,TMath::Power(10,xx), yy);
    gr->SetPointError(ip, TMath::Power(10, xx) - TMath::Power(10, xx-exl), TMath::Power(10, xx+exh)-TMath::Power(10, xx), eyl, eyh);
  }
}



Int_t SysUtils::GetCommonLimits(const Int_t ngr,  const TGraphAsymmErrors * grs[], Int_t id0s[], Int_t id1s[])
{
  Double_t maxlowerpoint= -1e2;
  Double_t minupperpoint = 1e2;
  for(Int_t ii=0; ii<ngr; ii++){
    Double_t gmin = -999;
    Double_t gmax = -999;
    Bool_t klogx = kFALSE;
    if(!grs[ii]->GetN()){
      printf("SysUtils::GetCommonLimits::bad graph n=0 igr %d %s %s\n", ii, grs[ii]->GetName(), grs[ii]->GetTitle());
      return 0;
      //exit(1);                                                                                                                                                                                    
    }
    GetGrPointLimits(grs[ii], gmin, gmax, klogx);
    if(gmin>maxlowerpoint){
      maxlowerpoint = gmin;
    }
    if(gmax<minupperpoint){
      minupperpoint = gmax;
    }
    printf("SysUtils::GetCommonLimits::loop in graph igr %d/%d min %e max %e maxlowerpoint %e minupperpoint %e\n", ii, ngr, gmin, gmax, maxlowerpoint, minupperpoint);
  }


  Int_t npoint = -999;
  for(Int_t igr=0; igr<ngr; igr++){
    id0s[igr]=-999;
    id1s[igr]=-999;
    for(Int_t pp=0; pp<grs[igr]->GetN(); pp++){
      const Double_t xx = grs[igr]->GetX()[pp];
      if(fabs(xx-maxlowerpoint)<EPSILON){
        id0s[igr] = pp;
      }
      if(fabs(xx-minupperpoint)<EPSILON){
        id1s[igr] = pp;
      }
    }
    if(id0s[igr]<0){
      printf("SysUtils::GetCommonLimits::not found for %s igr %d maxlowerpoint %e , %e %e\n", grs[igr]->GetName(), igr, maxlowerpoint, grs[igr]->GetX()[0], grs[igr]->GetX()[0]-maxlowerpoint); exit(1);
    }
    if(id1s[igr]<0){
      printf("SysUtils::GetCommonLimits::not found for igr %d minupperpoint %e\n", igr, minupperpoint);
      for(Int_t pp=0; pp<grs[igr]->GetN(); pp++){
        const Double_t xx = grs[igr]->GetX()[pp];
        printf("pp %d xx %e del %e\n", pp, xx, xx-minupperpoint);
      }
      exit(1);
    }
    printf("SysUtils::GetCommonLimits::results igr %d start %d maxlowerpoint %e end %d minupperpoint %e\n", igr, id0s[igr], maxlowerpoint, id1s[igr], minupperpoint);

    const Int_t tmpnp = id1s[igr]-id0s[igr]+1;
    if(npoint<0){
      npoint = tmpnp;
    }
    else if(tmpnp!=npoint){
      printf("SysUtils::GetCommonLimits::npoint not stable!! igr %d %d %d\n", igr, tmpnp, npoint);
      return -999;
      //exit(1);                                                                                                                                                                                    
    }
  }
  return npoint;
}


TGraphAsymmErrors *SysUtils::HistToGraph(const TString grname, const TH1 *hh, const Double_t thres, const TH1 *herr, const Double_t xmin, const Double_t xmax, const Bool_t klogx)
{
  const Int_t nbin = hh->GetNbinsX();
  Double_t xxs[nbin], yys[nbin], exls[nbin], exhs[nbin], eys[nbin];
  Int_t np=0;
  for(Int_t ii=1; ii<=nbin; ii++){
    const Double_t iyy = hh->GetBinContent(ii);
    if(iyy<=thres)
      continue;

    const Double_t iey = hh->GetBinError(ii);
    if(iey<EPSILON){
      if(iyy>EPSILON){
        printf("SysUtils::HistToGraph warning! should be fine if this is ratio %d %e %e\n", ii, iyy, iey); //exit(1);                                                                               
      } 
      //continue;                                                                                                                                                                                   
    }

    const Double_t ixx = klogx? TMath::Sqrt(hh->GetXaxis()->GetBinLowEdge(ii)*hh->GetXaxis()->GetBinUpEdge(ii)): hh->GetBinCenter(ii);
    if(ixx<xmin || ixx>xmax){
      //printf("test SysUtils::HistToGraph rejecting ixx %e xmin %e xmax %e\n", ixx, xmin, xmax);                                                                                                   
      continue;
    }

    Double_t iexl = 0;
    Double_t iexh = 0;
    if(herr){
      iexl = iexh = herr->GetBinContent(herr->GetXaxis()->FindBin(ixx));
    }
    else{
      if(klogx){
        iexl = ixx - hh->GetXaxis()->GetBinLowEdge(ii);
        iexh = hh->GetXaxis()->GetBinUpEdge(ii) - ixx;
      }
      else{
        iexl = iexh = hh->GetBinWidth(ii)/2.;
      }

    }

    xxs[np] = ixx;
    yys[np] = iyy;
    exls[np] = iexl;
    exhs[np] = iexh;
    eys[np] = iey;
    np++;
  }
  TGraphAsymmErrors * gr = new TGraphAsymmErrors(np, xxs, yys, exls, exhs, eys, eys);
  gr->SetName(grname);
  gr->SetMaximum(hh->GetMaximum());
  gr->SetMinimum(hh->GetMinimum());
  gr->GetXaxis()->SetLimits(hh->GetXaxis()->GetXmin(), hh->GetXaxis()->GetXmax());
  return gr;
}

TGraphAsymmErrors * SysUtils::TreeToGraph(const TString tname, const TString cut, Char_t * sx, Char_t * sy, Char_t * sex0, Char_t * sex1, Char_t * sey0, Char_t * sey1, const Bool_t kXLogToLinear, const Bool_t kprint)
{
  TTree * rawtree=(TTree*)gDirectory->Get(tname);
  if(!rawtree){
    printf("SysUtils::TreeToGraph no tree %s\n", tname.Data());exit(1);
  }

  TTree *cuttree = rawtree->CopyTree(cut);
  TTree * tmptree=(TTree*) cuttree->Clone(Form("tmptree%s%s",tname.Data(), sy));//necessary to avoid address confusion when executed repeatedly.                                                    

  Double_t vx = -1e10, vy = -1e10, vex0 =-1e10, vex1 = -1e10, vey0 = -1e10, vey1 = -1e10;
  tmptree->SetBranchAddress(sx, &vx);
  tmptree->SetBranchAddress(sy, &vy);
  if(strlen(sex0) && strlen(sex1)){
    tmptree->SetBranchAddress(sex0, &vex0);
    tmptree->SetBranchAddress(sex1, &vex1);
  }  
  if(strlen(sey0)){
    tmptree->SetBranchAddress(sey0, &vey0);
  }
  //if sey0 and sey1 are the same, vey0 will not be set from tree, has to be set by vey1                                                                                                            
  if(strlen(sey1)){
    tmptree->SetBranchAddress(sey1, &vey1);
  }

  TGraphAsymmErrors * tmpgr = new TGraphAsymmErrors;
  const Int_t ntot = tmptree->GetEntries();
  for(Int_t ii=0; ii<ntot; ii++){
    tmptree->GetEntry(ii);

    //remove trivial -1e10                                                                                                                                                                          
    if(vex0<-1e9) vex0=0;
    if(vex1<-1e9) vex1=0;
    if(vey0<-1e9) vey0=0;
    if(vey1<-1e9) vey1=0;

    //minos error low < 0 by convention                                                                                                                                                             
    vex0=TMath::Abs(vex0);
    vex1=TMath::Abs(vex1);
    vey0=TMath::Abs(vey0);
    vey1=TMath::Abs(vey1);

    if(strcmp(sey0,sey1)==0){
      vey0=vey1;
    }

    if(kXLogToLinear){
      const Double_t xtrue = pow(10, vx);
      const Double_t xlow = pow(10, vx-vex0);
      const Double_t xhigh = pow(10, vx+vex1);
      vx = xtrue;
      vex0 = xtrue - xlow;
      vex1 = xhigh - xtrue;
    }

    tmpgr->SetPoint(ii, vx, vy);
    tmpgr->SetPointError(ii, vex0, vex1, vey0, vey1);
    if(kprint){
      printf("SysUtils::TreeToGraph %s %s %s -- %d %e %e -- %e %e %e %e\n", tname.Data(), sx, sy, ii, vx, vy,  vex0, vex1, vey0, vey1);
    }
  }
  delete tmptree;

  TGraphAsymmErrors * gr = (TGraphAsymmErrors *)tmpgr->Clone(Form("gr%s%s",tname.Data(), sy));
  delete tmpgr;
  return gr;
}


void SysUtils::FitSlicesY(const TH2D *hh, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, TH1D *&hchi, const Double_t thres, const Double_t lowfrac, const Double_t highfrac, TList *ll)
{
  //                                                                                                                                                                                                
  //ffunc 3 parameter function: ALandau, AGaus                                                                                                                                                      
  //                                                                                                                                                                                                

  const Int_t x0 = hh->GetXaxis()->GetFirst();
  const Int_t x1 = hh->GetXaxis()->GetLast();
  const Int_t y0 = hh->GetYaxis()->GetFirst();
  const Int_t y1 = hh->GetYaxis()->GetLast();

  const Int_t nx = hh->GetNbinsX();
  const Int_t ny = hh->GetNbinsY();
  const Double_t xmin = hh->GetXaxis()->GetXmin();
  const Double_t xmax = hh->GetXaxis()->GetXmax();
  const Double_t ymin = hh->GetYaxis()->GetXmin();
  const Double_t ymax = hh->GetYaxis()->GetXmax();

  hnor = new TH1D(Form("%s_amp",hh->GetName()), "", nx, xmin, xmax);
  hmpv = new TH1D(Form("%s_mpv",hh->GetName()), "", nx, xmin, xmax);
  hwid = new TH1D(Form("%s_wid",hh->GetName()), "", nx, xmin, xmax);
  hres = new TH1D(Form("%s_res",hh->GetName()), "", nx, xmin, xmax);
  hchi = new TH1D(Form("%s_chi",hh->GetName()), "", nx, xmin, xmax);

  const Double_t *hxbins = GetAxisArray(hh->GetXaxis());

  hnor->GetXaxis()->Set(nx, hxbins);
  hmpv->GetXaxis()->Set(nx, hxbins);
  hwid->GetXaxis()->Set(nx, hxbins);
  hres->GetXaxis()->Set(nx, hxbins);
  hchi->GetXaxis()->Set(nx, hxbins);

  for(Int_t ix=x0; ix<=x1; ix++){
    TH1D *htmp = new TH1D(Form("%s_%d", hh->GetName(), ix),"",ny, ymin, ymax);
    //checked, ok                                                                                                                                                                                   
    const Double_t *hhybins = GetAxisArray(hh->GetYaxis());
    //hh->GetYaxis()->GetXbins()->GetArray();                                                                                                                                                       
    if(hhybins){
      htmp->GetXaxis()->Set(ny, hhybins);
    }

    Double_t ntot = 0;
    for(Int_t iy=y0; iy<=y1; iy++){
      const Double_t be = hh->GetBinError(ix,iy);
      const Double_t bc = hh->GetBinContent(ix, iy);

      if(be<EPSILON){
        if(bc>EPSILON){
          printf("SysUtils::FitSlicesY error %d %d %e %e\n", ix, iy, be, bc); exit(1);
        }
        continue;
      }

      htmp->SetBinContent(iy, bc);
      htmp->SetBinError(iy, be);

      ntot += (bc/be)*(bc/be);

      //if(be) printf("test %d %d : %f %f %f\n", ix, iy, bc, be, pow(bc/be,2));                                                                                                                     
    }

    hnor->SetBinContent(ix, ntot);
    hnor->SetBinError(  ix, 0);

    if(ntot<thres || htmp->GetRMS()<EPSILON){
      delete htmp;
      continue;
    }

    //test htmp->Draw();                                                                                                                                                                            
    Double_t pars[10]={htmp->Integral(0,htmp->GetNbinsX()+1)*htmp->GetBinWidth(1), htmp->GetMean(), htmp->GetRMS()};
    Double_t errs[10], chi[10];

    TH1D * tmpslicefit = 0x0;
    pars[1] = TruncatedMean(htmp, lowfrac, highfrac, &(pars[2]), &(errs[1]));
    errs[2]=0;
    chi[0] = chi[1] = 0;

    pars[2]=TMath::Abs(pars[2]);
    //hnor->SetBinContent(ix, htmp->GetBinContent(htmp->GetMaximumBin()));//htmp->Integral(0,htmp->GetNbinsX()+1));                                                                                 
    hmpv->SetBinContent(ix, pars[1]);
    hmpv->SetBinError(  ix, errs[1]);

    hwid->SetBinContent(ix, pars[2]);
    hwid->SetBinError(  ix, errs[2]);

    hres->SetBinContent(ix, fabs(pars[1])>EPSILON? pars[2]/fabs(pars[1]):0);
    hres->SetBinError(  ix, fabs(pars[1])>EPSILON? errs[2]/fabs(pars[1]):0);

    hchi->SetBinContent(ix, chi[1]>=1 ? chi[0]/chi[1]: 0);
    hchi->SetBinError(ix, 0);

    if(ll){
      ll->Add(htmp);
      if(tmpslicefit)
        ll->Add(tmpslicefit);
    }
    else{
      delete htmp;
      delete tmpslicefit;
    }
  }

  TH1 *hhs[]={hnor, hmpv, hwid, hres, hchi};
  const TString yt[]={"N", "MPV", "#sigma", "#sigma/MPV", "#chi^{2}/NDOF"};
  const Int_t nh = sizeof(hhs)/sizeof(TH1*);
  for(Int_t ii=0; ii<nh; ii++){
    hhs[ii]->SetYTitle(Form("%s of %s", yt[ii].Data(), hh->GetYaxis()->GetTitle()));
    hhs[ii]->SetXTitle(hh->GetXaxis()->GetTitle());
    hhs[ii]->GetYaxis()->SetTitleOffset(hh->GetYaxis()->GetTitleOffset());
    hhs[ii]->SetTitle(hh->GetTitle());
  }
}


Double_t SysUtils::ToPt(const Double_t mom, const Double_t eta)
{
  /*                                                                                                                                                                                                
    TLorentzVector b                                                                                                                                                                                
    b.SetPtEtaPhiE(0.567,0.123,0,0)                                                                                                                                                                 
    b.Pt()-SysUtils::ToPt(b.P(), b.Eta())                                                                                                                                                           
    //should be eps                                                                                                                                                                                 
    */

  const Double_t sintheta = TMath::Sin( 2*TMath::ATan( TMath::Exp( TMath::Abs( eta ) ) ) );
  return mom*sintheta;
}


void SysUtils::ScaleGraphPoint(TGraphAsymmErrors * gr, const Int_t ipoint, const Double_t fac)
{
  const Double_t xx = gr->GetX()[ipoint];
  const Double_t yy = gr->GetY()[ipoint];
  const Double_t el = gr->GetEYlow()[ipoint];
  const Double_t eh = gr->GetEYhigh()[ipoint];
  gr->SetPoint(ipoint,xx, yy*fac);
  gr->SetPointEYlow(ipoint, el*fac);
  gr->SetPointEYhigh(ipoint, eh*fac);
}


TGraphAsymmErrors * SysUtils::MeasurementConversion(const TString hn, const TGraphAsymmErrors * gv1, const TH2D * hy2x1, const Bool_t kcutEdge, const Double_t thres)
{
  
  const Int_t matrixFirstX = hy2x1->GetXaxis()->GetFirst();
  //check nbins

  const Int_t nv1 = gv1->GetN();
  const Int_t matrixNx = hy2x1->GetXaxis()->GetLast() - matrixFirstX + 1;
  if(matrixNx!=nv1){
    printf("SysUtils::MeasurementConversion hy2x1 nv1 != hy2xinx %s: %d , %s %d\n", gv1->GetName(), nv1, hy2x1->GetName(), matrixNx); return 0x0;
  }
  
  //check v1 binning

  for(Int_t ii=0; ii<nv1; ii++){
    if( fabs( gv1->GetX()[ii] - hy2x1->GetXaxis()->GetBinCenter(ii+matrixFirstX) ) > EPSILON ){
      printf("SysUtils::MeasurementConversion gv1 hy2x1 binning does NOT match at ii %d gv1 %e hy2x1 %e\n", ii, gv1->GetX()[ii], hy2x1->GetXaxis()->GetBinCenter(ii+matrixFirstX)); return 0x0;
    }
  }
  
  const Int_t nv2 = hy2x1->GetNbinsY();
  TGraphAsymmErrors * gout = new TGraphAsymmErrors;
  gout->SetName(hn);
  gout->SetTitle(Form("converted from %s %s", gv1->GetName(), hy2x1->GetName()));
  
  for(Int_t i2 = 0; i2< nv2; i2++){
    Double_t ff = 0;
    Double_t sl2 = 0;
    Double_t sh2 = 0;

    //Double_t Integral(Int_t binx1, Int_t binx2, Int_t biny1, Int_t biny2, Option_t* option = "") const

    const Double_t nor = hy2x1->Integral(matrixFirstX, nv1+matrixFirstX-1, i2+1, i2+1);

    //if on the edge there is counts, it means the acceptance is not complete

    const Double_t edgeThres = nor*1e-2;

    //the bin is large enough, no need to have two bins

    //const Double_t edgeCount = hy2x1->Integral(matrixFirstX, matrixFirstX+1, i2+1, i2+1) + hy2x1->Integral(nv1+matrixFirstX-2, nv1+matrixFirstX-1, i2+1, i2+1);

    const Double_t edgeCount = hy2x1->GetBinContent(matrixFirstX,i2+1) + hy2x1->GetBinContent(nv1+matrixFirstX-1, i2+1);


    if(nor>thres){
      if(!kcutEdge || edgeCount< edgeThres){

        for(Int_t i1 = 0; i1 < nv1; i1++){
          const Double_t v1f = gv1->GetY()[i1];
          const Double_t ww = hy2x1->GetBinContent(i1+matrixFirstX, i2+1)/nor;
          
          if(v1f<EPSILON || ww<EPSILON){
            continue;
          }
          
          const Double_t v1el = gv1->GetEYlow()[i1];
          const Double_t v1eh = gv1->GetEYhigh()[i1];

          const Double_t dw = TMath::Sqrt(ww*(1-ww)/nor);

          //v1elh correspond to the whole pbin ~1/sqrt(totPstat)

          //(p,z) corresponds to only bin content ~1/sqrt(bincontent)

          const Double_t totPstat = hy2x1->Integral(i1+matrixFirstX, i1+matrixFirstX, 0, nv2);
          const Double_t inflat = totPstat/hy2x1->GetBinContent(i1+matrixFirstX, i2+1);

          ff += v1f*ww;
          //assuming uncorrelated errors

          sl2 += TMath::Power(v1f*ww, 2) * ( TMath::Power(dw/ww,2) + TMath::Power(v1el/v1f,2)*inflat );//ww*ww*v1el*v1el; ww*ww*v1eh*v1eh;
          sh2 += TMath::Power(v1f*ww, 2) * ( TMath::Power(dw/ww,2) + TMath::Power(v1eh/v1f,2)*inflat );//ww*ww*v1el*v1el; ww*ww*v1eh*v1eh;

          //printf("SysUtils::MeasurementConversion test %s i2 %d/%d yy %e nor %e , i1 %d/%d, v1f %e v1el %e v1eh %e ww*nor %e ff %e sl2 %e sh2 %e\n", gout->GetName(), i2, nv2, hy2x1->GetYaxis()->GetBinCenter(i2+1),  nor, i1, nv1, v1f, v1el, v1eh, ww*n\or, ff, sl2, sh2);
        }
        
        const Double_t efl = TMath::Sqrt(sl2);
        const Double_t efh = TMath::Sqrt(sh2);
        
        const Int_t np = gout->GetN();
        gout->SetPoint(np, hy2x1->GetYaxis()->GetBinCenter(i2+1), ff);
        gout->SetPointError(np, hy2x1->GetYaxis()->GetBinWidth(i2+1)/2, hy2x1->GetYaxis()->GetBinWidth(i2+1)/2, efl, efh);
        
      }
      else{
        printf("SysUtils::MeasurementConversion not calculated %s i2 %d nor %e kcutEdge %d edgeCount %e %e %e edgeThres %e\n", gout->GetName(), i2, nor, kcutEdge, hy2x1->GetBinContent(matrixFirstX,i2+1), hy2x1->GetBinContent(nv1+matrixFirstX-1, i2+1), edgeCount, edgeThres);
        //do not record any more if reaching edge
        if(gout->GetN()){
          break;
        }
      }
    }
  }
  
  return gout;
}


Bool_t SysUtils::GetGrPointLimits(const TGraphAsymmErrors * gr, Double_t & gmin, Double_t & gmax, Bool_t & klogx)
{
  if(!gr){
    printf("SysUtils::GetGrPointLimits gr null\n"); exit(1);
  }
  if(gr->GetN()==0){
    printf("SysUtils::GetGrPointLimits gr n=0!!\n"); exit(1);
  }

  gmin = 1e10;
  gmax = -1e10;
  for(Int_t ii=0; ii<gr->GetN(); ii++){
    const Double_t xx = gr->GetX()[ii];
    if(xx>gmax)
      gmax = xx;
    if(xx<gmin)
      gmin = xx;
  }

  const Double_t e0 = gr->GetEXlow()[0];
  const Double_t e1 = gr->GetEXhigh()[0];

  if(fabs(e0-e1)>EPSILON)
    klogx = kTRUE;
  else
    klogx = kFALSE;

  Bool_t kup=kTRUE;
  if( fabs(gr->GetX()[0]-gmin)>EPSILON && fabs(gr->GetX()[gr->GetN()-1]-gmax)<EPSILON){
    kup = kTRUE;
  }
  else{
    kup = kFALSE;
  }

  return kup;
}


Double_t * SysUtils::GetAxisArray(TAxis * aa)
{
  TArrayD * xs = new TArrayD(*(aa->GetXbins()));
  if(xs->GetSize()==0){
    const Int_t nbin=aa->GetNbins();
    xs->Set(nbin+1);
    for(Int_t ii=0; ii<=nbin; ii++){
      xs->SetAt(aa->GetBinUpEdge(ii),ii);
    }
  }

  Double_t * bins = xs->GetArray();
  if(!bins){
    printf("SysUtils::GetAxisArray bins null!! %d\n", xs->GetSize());exit(1);
  }
  return bins;
}


Double_t SysUtils::GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms, Double_t * gerr)
{
  //                                                                                                                                                                                                                                                         
  //calculate mean (with error) and rms from sum, w2s, nn                                                                                                                                                                                                    
  //if nn=0, mean, error, and rms all = 0                                                                                                                                                                                                                    
  //                                                                                                                                                                                                                                                         

  Double_t tmean = 0, trms = 0, terr = 0;

  const Double_t eps = 1e-10;
  if(nn>eps){
    tmean = sum/nn;

    const Double_t arg = w2s/nn-tmean*tmean;
    if(TMath::Abs(arg)<eps){
      trms = 0;
    }
    else{
      if( arg <0 ){
        printf("SysUtils::GetMeanRMS error negative sqrt argument!! %e -- %e %e %f\n", arg, w2s, sum, nn); exit(1);
      }

      trms = TMath::Sqrt(arg);
    }

    terr = trms/TMath::Sqrt(nn);
  }

  if(grms){
    (*grms) = trms;
  }

  if(gerr){
    (*gerr) = terr;
  }

  return tmean;

}

Double_t SysUtils::TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr)
{
  //                                                                                                                                                                                                                                                         
  //do truncation on histogram                                                                                                                                                                                                                               
  //                                                                                                                                                                                                                                                         
  //if hh is scaled, be sure Sumw2 is called before scaling!! then mean, rms and err will all be correct                                                                                                                                                     

  //with under-/over-flow                                                                                                                                                                                                                                    
  Double_t npreTrunc = 0;
  for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
    const Double_t be = hh->GetBinError(itmp);
    const Double_t bc = hh->GetBinContent(itmp);
    if(be<EPSILON){
      if(bc>EPSILON){
        printf("SysUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
      }
      continue;
    }
    npreTrunc += bc*bc/be/be;
  }

  const Double_t nstart = npreTrunc*lowfrac;
  const Double_t nstop = npreTrunc*highfrac;

  //with Double_t this should also handle normalized hist
                                                                                                                                                                                                    
  Double_t ntot = 0;
  Double_t nused = 0;
  Double_t sum = 0;
  Double_t w2s = 0;
  for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
    const Double_t be = hh->GetBinError(itmp);
    const Double_t bc = hh->GetBinContent(itmp);
    if(be<EPSILON){
      if(bc>EPSILON){
        printf("SysUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
      }
      continue;
    }
    const Double_t weight = bc*bc/be/be;
    ntot+=weight;
    //<= correct, because when high=1, nstop has to be included 

    if(ntot>nstart && ntot<=nstop){

      const Double_t val = hh->GetBinCenter(itmp);
      sum += weight*val;
      w2s += weight*val*val;

      nused += weight;

      //printf("test %d %f %f --- %f %f -- %f %f\n", itmp, weight, val, sum, w2s, nused, nsample); 
    }
    else if(ntot>nstop){
      if(itmp>=hh->GetNbinsX()){
        printf("SysUtils::TruncatedMean warning hist range too small %s %f %f %d %d, %15f %15f %15f; nused w2s sum set to 0\n", hh->GetName(), hh->GetBinLowEdge(1), hh->GetBinLowEdge(itmp), itmp, hh->GetNbinsX(), hh->GetBinContent(hh->GetNbinsX())/hh->Integral(0,hh->GetNbinsX()+1), hh->GetBinContent(hh->GetNbinsX()), hh->Integral(0,hh->GetNbinsX()+1)); //exit(1);

        nused = 0;
        w2s = sum = 0;
      }
      break;
    }
  }

  return GetMeanRMS(nused, sum, w2s, grms, gerr);

}


void SysUtils::GetGrSysMode0(TList *ll, const TString tag, const Bool_t *badgr, const Bool_t kfit, const Double_t hmin, const Double_t hmax)
{
  printf("\n\n********************** SysUtils::GetGrSysMode0 ********************** \n\n");
  ll->ls();

  const Int_t rawn = ll->GetEntries();

  const TGraphAsymmErrors * grall[rawn];
  Int_t id0s[rawn];
  Int_t id1s[rawn];
  for(Int_t ii=0; ii<rawn; ii++){
    grall[ii] = 0x0;
    id0s[ii]=id1s[ii]=-999;
  }
  Int_t ng = 0;
  for(Int_t ii=0; ii<rawn; ii++){
    if(badgr && badgr[ii]){
      //skip bad gr
      printf("SysUtils::GetGrSysMode0 %s gr %d finally-kicked-out!\n", tag.Data(), ii);
      continue;
    }
    grall[ng]=(TGraphAsymmErrors*) ll->At(ii);
    ng++;
  }

  const Int_t npoint = GetCommonLimits(ng, grall, id0s, id1s);
  if(npoint<0){
    printf("SysUtils::GetGrSysMode0 npoint<0 %d GetCommonLimits not stable?\n", npoint); exit(1);
  }
  printf("SysUtils::GetGrSysMode0 ng %d np %d\n", ng, npoint);

  //output to ll by adding in the very end, not to disturb ng 

  TGraphAsymmErrors * grstat = new TGraphAsymmErrors;
  grstat->SetName(Form("grstat%s", tag.Data()));

  TGraphAsymmErrors * grsys = new TGraphAsymmErrors;
  grsys->SetName(Form("grsys%s", tag.Data()));
  grsys->SetFillStyle(3004);
  grsys->SetFillColor(kGray);

  ll->Add(grstat);
  ll->Add(grsys);

  if(!npoint){
    return;
  }

  //==================================================

  Double_t xx[npoint];
  Double_t exl[npoint];
  Double_t exh[npoint];
  Double_t yy[npoint];
  Double_t sumw2y[npoint];

  Double_t eyl[npoint];
  Double_t eyh[npoint];
  Double_t ymin[npoint];
  Double_t ymax[npoint];
  Int_t gmin[npoint];
  Int_t gmax[npoint];
  TH1D *hdists[npoint];

  TH2D * h2d=new TH2D(Form("h2d%s",tag.Data()),"",50, 0.1, 100, 100, hmin, hmax); ll->Add(h2d);
  BinLogX(h2d->GetXaxis());

  Int_t nread = 0;

  const Int_t hnbx = 2000;

  for(Int_t ix=0; ix<npoint; ix++){
    xx[ix]=-999;
    exl[ix]=0;
    exh[ix]=0;
    yy[ix]=0;
    sumw2y[ix]=0;
    eyl[ix]=0;

    ymin[ix]= 1e10;
    ymax[ix]=-1e10;
    gmin[ix]=-1;
    gmax[ix]=-1;

    hdists[ix] = 0x0;

    //=============================
    //loop over graphs

    Double_t grstaterrl[ng];
    Double_t grstaterrh[ng];

    for(Int_t ig=0; ig<ng; ig++){
      grstaterrl[ig] = 1E10;
      grstaterrh[ig] = 1E10;

      const TGraphAsymmErrors * gr = grall[ig];

      if(gr->GetN()<npoint){
        printf("SysUtils::GetGrSysMode0 ig %d/%d gr->GetN()<npoint %d %d\n", ig, ng, gr->GetN(), npoint); exit(1);
      }

      const Int_t iraw = id0s[ig]+ix;

      Double_t tmpx=-999, tmpy=-999;
      gr->GetPoint(iraw, tmpx, tmpy);
      printf("SysUtils::GetGrSysMode0 reading %s %d %e %e\n", gr->GetName(), iraw, tmpx, tmpy);

      if(xx[ix]<-990){
        xx[ix]=tmpx;
        exl[ix]=gr->GetEXlow()[iraw];
        exh[ix]=gr->GetEXhigh()[iraw];
      }
      else{
        const Double_t eps = EPSILON;
        if(TMath::Abs(xx[ix]-tmpx)>eps){
          printf("SysUtils::GetGrSysMode0 xx[%d]!=tmpx %d, %e %e %e %e\n", ix, iraw, xx[ix], tmpx, TMath::Abs(xx[ix]-tmpx), eps);exit(1);
        }
      }

      if(!hdists[ix]){
        hdists[ix] = new TH1D(Form("hdist%s_%d", tag.Data(), ix), Form("xx %e", tmpx), hnbx, hmin, hmax);    ll->Add(hdists[ix]);
      }

      h2d->Fill(tmpx, tmpy);
      hdists[ix]->Fill(tmpy);

      yy[ix]+=tmpy;
      sumw2y[ix]+=tmpy*tmpy;

      //eyl[ix]+= gr->GetEYlow()[iraw]; 
      //eyh[ix]+= gr->GetEYhigh()[iraw]; 

      grstaterrl[ig] = gr->GetEYlow()[iraw];
      grstaterrh[ig] = gr->GetEYhigh()[iraw];

      if(tmpy>ymax[ix]){
        ymax[ix]=tmpy;
        gmax[ix]=ig;
      }

      if(tmpy<ymin[ix]){
        ymin[ix]=tmpy;
        gmin[ix]=ig;
      }

      if(ix==0){
        nread++;
      }
    }
    //Double_t SysUtils::TruncatedMean(const Int_t nx, const Double_t xdata[], const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr, Double_t *wws) 
    //printf("testtest ng %d \n", ng); 
                                                                                                                                                                                                       
    for(Int_t itest=0; itest<ng; itest++){
      if(TMath::IsNaN(grstaterrl[itest])) grstaterrl[itest]=1e3;
      if(TMath::IsNaN(grstaterrh[itest])) grstaterrh[itest]=1e3;

      //printf("testtest %d %e %e\n", itest, grstaterrl[itest], grstaterrh[itest]); 
    }
    const Double_t lowf = 0.05;
    const Double_t highf = 0.2;

    const Int_t itmpstart = (Int_t)(lowf*ng);
    const Int_t itmpstop  = (Int_t)(highf*ng);
    if(itmpstop-itmpstart>0){
      eyl[ix] = TruncatedMean(ng, grstaterrl, lowf, highf);
      eyh[ix] = TruncatedMean(ng, grstaterrh, lowf, highf);
    }
    else{
      if(ng!=1){
        printf("SysUtils::GetGrSysMode0 wrong ng lowf highf %d %e %e\n", ng, lowf, highf);exit(1);
      }
      else{
        eyl[ix]= grstaterrl[0];
        eyh[ix]= grstaterrh[0];
      }
    }
  }  
  for(Int_t ii=0; ii<npoint; ii++){

    Double_t pars[10]={ hdists[ii]->Integral(0,hdists[ii]->GetNbinsX()+1)*hdists[ii]->GetBinWidth(1), hdists[ii]->GetMean(), hdists[ii]->GetRMS() };

    if(kfit){
    }

    const Double_t fy = pars[1];
    const Double_t rmsy = TMath::Abs(pars[2]);

    const Double_t sysel= rmsy; //fy-ymin[ii];
    const Double_t syseh= rmsy; //ymax[ii]-fy;

    //they are already the truncated mean 0.1-0.9  

    const Double_t statel = eyl[ii]; //eyl[ii]/nread; 
    const Double_t stateh = eyh[ii]; //eyh[ii]/nread; 

    printf("SysUtils::GetGrSysMode0 setting %d %e %e %e %e -- %e %e -- gmin %d gmax %d\n", ii, xx[ii], fy, sysel, syseh, ymin[ii], ymax[ii], gmin[ii], gmax[ii]);

    grstat->SetPoint(ii, xx[ii], fy);
    grstat->SetPointError(ii, exl[ii], exh[ii], statel, stateh);

    grsys->SetPoint(ii, xx[ii], fy);
    grsys->SetPointError(ii, exl[ii], exh[ii], sysel, syseh);
  }
}


void SysUtils::BinLogX(TAxis *axis)
{
  //                                                                                                                                                                                                                                          
  // Method for the correct logarithmic binning of histograms                                                                                                                                                                                 
  // copied and modified from AliTPCcalibBase                                                                                                                                                                                                 

  const Int_t bins = axis->GetNbins();

  const Double_t from = axis->GetXmin();
  const Double_t to = axis->GetXmax();
  if (from<EPSILON) return;
  Double_t *new_bins = new Double_t[bins + 1];

  new_bins[0] = from;
  const Double_t factor = TMath::Power(to/from, 1./bins);

  for (int i = 1; i <= bins; i++) {
    new_bins[i] = factor * new_bins[i-1];
  }
  axis->Set(bins, new_bins);
  delete [] new_bins;
}



Double_t SysUtils::TruncatedMean(const Int_t nx, const Double_t xdata[], const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr, Double_t *wws)
{

  Int_t index[nx];
  TMath::Sort(nx, xdata, index, kFALSE);

  Int_t nused = 0;
  Double_t sum = 0;
  Double_t w2s = 0;
  const Int_t istart = Int_t (nx*lowfrac);
  const Int_t istop  = Int_t (nx*highfrac);

  //=,< correct, because when low=0, high=1 it is correct                                                                                                                                                                                     
  for(Int_t ii=istart; ii<istop; ii++){
    Double_t weight = 1;
    if(wws){
      weight = wws[index[ii]];
    }
    const Double_t sx = xdata[index[ii]]*weight;

    sum += sx;
    w2s += sx*sx;

    nused++;
    //printf("test in loop %d/%d %f %f %f\n", ii, nused, sx, sum, w2s);                                                                                                                                                                       

  }

  return GetMeanRMS(nused, sum, w2s, grms, gerr);

}

#endif


